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We introduce a unitary matrix-product operator (UMPO) based variational method that approximately finds 
all the eigenstates of fully many-body localized (fMBL) one-dimensional Hamiltonians. The computational cost 
of the variational optimization scales linearly with system size for a fixed bond dimension of the UMPO ansatz. 
We demonstrate the usefulness of our approach by considering the Heisenberg chain in a strongly disordered 
magnetic field for which we compare the approximation to exact diagonalization results. 


Introduction: The phenomenon of many-body localization 
(MBL) generalizes Anderson localization m (AL) to interact¬ 
ing systems EM). In the Anderson problem the many-body 
Fock/SIater states constructed from the single particle states 
have two important features. First, they exhibit an economi¬ 
cal description —L single particle states for a system of size 
L are sufficient to construct all 2^ many-body states. Second, 
all many-body states exhibit an area law for the entanglement 
entropy stemming from the localized nature of the constituent 
single particle states. Naturally, attention has focused on what 
happens to these two properties in the MBL regime. 

It was noted early on iQ that many-body eigenstates in the 
MBL regime would have only local entanglement and thus 
obey area laws iHS). Subsequently Bauer and Nayak m ex¬ 
amined the behavior of the entanglement entropy in detail and 
demonstrated the area law as well as deviations due to rare 
regions and states. In another set of papers MM the phe¬ 
nomenology of MBL systems was traced to an emergent set of 
L commuting local integrals of motion (often called “I-bits”) 
which are believed to exist in fMBL systems—i.e. systems in 
which all many-body eigenstates are localized. 

These two developments invite a natural closure in which 
the full set of 2^ many-body eigenstates are explicitly con¬ 
structed from 0{L) local ingredients, at least approximately. 
The well known connection of the area law to matrix-product 
state (MPS)/tensor network representations of many-body 
states suggests that the latter are the correct language in which 
to carry out this program. The program has two components: 
showing that such a compact representation exists and pro¬ 
viding a recipe for finding it without recourse to a knowledge 
of the exact eigenstates, potentially rendering a much larger 
range of system sizes computationally tractable. 

In an important development, two groups have addressed 
the existence problem. Building on earlier work ca, Pekker 
and Clark (PC) Ca have shown that the unitary operators that 
exactly diagonalize fMBL systems can be represented by ma¬ 
trix products operators (MPOs) ifT^ of bond dimensions that 
appear to grow very slowly with system size—in contrast to 
delocalized systems where the dimension grows exponentially 
with system size. The slow growth that they do observe is pre¬ 
sumably due to rare many-body resonances/Griffiths effects; 
in its absence, the MPOs would yield the sought after 0(L) lo- 
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FIG. I. (a) Schematic representation of an MPS representation of a 
state 1^). (b) Variational ansatz for the unitary U that encodes all 
eigenstates of a fully many-body localized Hamiltonian. The local 

unitaries are parametrized as with real symmet¬ 

ric matrices , n = I... L — 1 and m = I... A^iayer. 


cal description of the full spectrum. Parallel work ca argued 
for the congruent result that the presence of local integrals of 
motion implies the existence of a single “spectral tensor net¬ 
work” that efficiently represents the entire spectrum of energy 
eigenstates in the fMBL phase. These developments however 
have not led to an algorithm for finding a compact represen¬ 
tation directly and even finding MPOs representing exactly 
known diagonalizing unitaries a la PC scales exponentially 
with system size CS). 

In this paper we propose an approach to directly and ef¬ 
ficiently find an approximate compact representation of the 
diagonalizing unitary by using a variational unitary MPO 
(VUMPO) ansatz. To this end, we construct a cost function 
whose minimum yields the exact unitary and, hence, the en¬ 
tire set of 2^ exact eigenstates of a system of L qubits. We 
show that for a fixed bond dimension of the approximate U, 
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optimizing the cost-function in d = 1 can be performed at a 
computational cost that is only linear in system size which, 
in theory, allows us to access system sizes far beyond those 
possible by ED. 

MPS and MPO notation: An MPS representation of a quan¬ 
tum state living in a basis spanned by L qubits takes the form 

i^) = E E (1) 

{a-} 0<ji<D 

whereas an MPO representation of an operator in the same 
Hilbert space takes the form 

( 2 ) 

where cr^, G {t, i} and we use a compact notation in which 
a = cTi, (72, • • • , ctl denotes the 2^ states (analogous for 
r). Figureshows a pictorial representation of these objects. 
The MPSs/MPOs are represented by rank three/four tensors 
on each site i (except the first and last tensors which 
are rank two/three); the external leg(s) (J^, refer to the phys¬ 
ical spin indices whereas the ji are the internal virtual indices 
that are contracted. Each is a D‘^ dimensional 

matrix where D is the bond-dimension of the matrix. 
Method: We now introduce the VUMPO ansatz and an algo¬ 
rithm to numerially optimize it. Let us assume that H is an 
fMBL Hamiltonian defined on an L-site chain of spin 1/2 op¬ 
erators. It is our goal to find a unitary MPO approximation U 
of the unitary that diagonalizes the Hamiltonian such that the 
2^ eigenstates of H are given by 

(3) 

In the parlance of Refs. (ToKTTl, the physical basis operators 
ai are the “p-bits” wheras the are the local “1-bits”. Each 
eigenstate is labeled by the occupation of l-bits r = {tti 
• • • t}, and is obtained by acting with the MPO representation 
of U on the product state |t). In this language of MPOs, it is 
clear how the 2^ MB eigenstates are constructed from the L 
matrices ; further, if the bond-dimension of the matrices 
scales as 0(1) with the system size, the eigenstates are only 
locally entangled in the p-bit basis and a description of the full 
eigenbasis in terms of 0{L) local ingredients is possible. 

The VUMPO is found by numerically minimizing the cost 
functional 

/({AW}) = > 0, (4) 

with The cost function is the vari¬ 

ance of the energy summed over all approximate MB eigen¬ 
states. Naively, one might expect the time to evaluate the 
cost function Eq. 0 to scale exponentially with the sys¬ 
tem size L as the sum is performed over 2^ MB eigen¬ 
states. However, remarkably, the computation can be per¬ 
formed in a time scaling linearly with system size ifTHl ! 
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FIG. 2. Comparison of the exact energy levels (blue lines) with the 
ones found by the variational optimization (red lines) for VF = 8 and 
L = 8 as a function of the number of layers of two-site gates. The 
right panel shows a zoom of some energy levels at the bottom and in 
the center of the spectrum. 

For example, the term evaluated 

by “doubling” the degrees of freedom and defining a state 
\(j)) = With this notation we find that 

~ (^ 1 ^ ( 8 ) ( 8 ) 1|0). This expectation 

value can be efficiently evaluated using the MPO formalism 
and the most expensive part of the evaluation scales, for a 
given Hamiltonian in MPO form, as oc LD^ (see Appendix 
A for details and a diagrammatic representation of the terms). 
One can now iteratively minimize / by locally optimizing 
each using, for example, the conjugate gradient algo¬ 
rithm. 

In general, an MPO compression of a unitary operator will 
not strictly respect unitarity. To get a valid positive-definite 
cost function in these cases, we need to add a Lagrange mul¬ 
tiplier to enforce unitarity (or consider other cost functions 
which don’t assume orthonormality of the eigenstates). In 
practice, these methods lead to either very unstable, or very 
computationally expensive optimizations. 

The key to a stable optimization protocol turns on restrict¬ 
ing our algorithm to the manifold of strictly unitary MPOs of 
a given bond-dimension. To achieve this, we parameterize the 
VUMPO as a finite depth circuit of two-site unitaries as shown 
in Fig.[2b). This Ansatz incorporates two important proper¬ 
ties: (i) The VUMPO is unitary for all parameters and (ii) it 
is local for any finite Mayer- We use a single unitary to obtain 
all eigenstates, but readers will note the obvious connection to 
the quantum computational notion ITTl that each MBL eigen¬ 
state can be constructed from a reference Fock state via the 
operation of a, in general different, finite depth circuit made 
up of local unitaries. Finally, we note that we can rewrite the 
unitary network as a strictly unitary MPO with bond dimen¬ 
sion D < , where Mayer is the number of layers of 

two-site gatestm. However, this step is not necessary and we 
can evaluate the cost function by directly contracting the uni- 
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taries circuit which, in fact, gives a considerable speed up for 
the systems we consider here 113. 

The algorithm to find the VUMPO is then similar in spirit 
to the density matrix renormalization group (DMRG) method 
1^ , except instead of finding the lowest energy state, we min¬ 
imize the cost function Eq. Q by sweeping through the local 
unitaries: 

[ml 

(i) Initialize the local unitaries = e t^^i by choosing 

random symmetric matrices where n = 1, 2, • • • L and 

TTL 1, 2, • • • Allayer* 

(ii) Locally minimize the cost function by varying the ele¬ 
ments of a given by using, e.g., a conjugate gradient 
method. 

(iii) Update the network and repeat the previous step for the 
next unitary. 

(iv) Continue the sweeping procedure by minimizing the local 
unitaries successively until convergence. A full sweep across 
all the unitaries has to scale as 0{L). 

We find that the number of steps needed for convergence ap¬ 
pears to be approximately independent of L. This gives an 
overall scaling of the algorithm as 0{LD^) ^ 

Once the algorithm has converged, the VUMPO can be used 
to obtain all the eigenstates of the system, and to efficiently 
compute observables using the MPS formalism. 

Results: We consider the Heisenberg model with random 2;- 
directed magnetic fields: 

H = J^Sn-Sn+i-Y^ KSl (5) 

n n 

where Sn are spin 1/2 operators and the fields hn are drawn 
randomly from the interval [—lU, lU] and we set J = 1. This 
model has been studied extensively in the context of MBL and 
several numerical studies strongly suggest that H is fMBL for 
VE > 3.5 EHmEl. 

Energy Spectrum: We begin by comparing the energies ob¬ 
tained using the VUMPO approach with the exact spectrum 
(full diagonalization). The converged results for lU = 8 and 
L = S with different numbers of layers Mayer are shown in 
Fig. For Mayer = 0, the VUMPO is the identity (i.e, no 
variational parameters) and the resulting approximate eigen¬ 
states are simple product states of the form |cri)|cr 2 )... Ictl) 
with (In =t, The overall bandwidth in this case agrees rel¬ 
atively well with the exact results because W is the domi¬ 
nant energy scale in the problem. However, as shown in the 
zoomed in plots, the deviation of individual energy levels is 
relatively large compared to the mean-level spacing because 
the product states completely neglect local quantum fluctu¬ 
ations which are present in the exact eigenstates. Increasing 
Mayer Strongly improves the agreement between the exact and 
approximate energy levels as the network successively adds 
entanglement over longer distances. 

Next we turn to the mean variance of the energy, which 
is simply the disorder averaged cost function Eq. 0 divided 
by 2^. Figure shows this quantity disorder averaged over 
50 realizations as a function of system size for different fixed 



L 

FIG. 3. Mean variance of the energy as a function of system size for 
different number of layers for W — 8. 
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FIG. 4. Comparison of the exact spectral function A{cjj) (black dots) 
with those obtained using different approximations (see text for de¬ 
tails) for L = 10 and FF = 8. Spectra are shown using a Lorentzian 
broadening with an imaginary part of e = 0.1. Inset: Same data with 
W = 16. 


Mayer- Wc obscrvc a linear increase of the mean variance 
with system size, and find that the slope decreases as Mayer 
is increased. This tells us that for a given Mayer our approxi¬ 
mate eigenstates involve a constant error per unit length which 
decreases as Mayer is increased. This scaling is entirely ex¬ 
actly the same as that obtained for ground states obtained via 
DMRG. 

Spectral Functions: To examine the quality of our approx¬ 
imated eigenstates (with a view to capturing local properties), 
we use the VUMPO ansatz to obtain the infinite-temperature 
spectral function 


A{lo) 


^ E \{Tl\Sl/,\T2)\^6{u-E^,+Er,). (6) 


Spectral functions can again be efficiently evaluated using 
matrix-product techniques and it is also possible to effi¬ 
ciently target different energy densities by considering finite- 
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temperature spectral functions (Hi [23. Figure compares 
the A{uj) obtained using the VUMPO approach for L = 10 
with different disorder strengths and Mayer = 0,1, 2 with the 
exact results. The spectral functions are dominated by a large 
peak at cc = 0 which reflects the strongly localized nature of 
the eigenstates, i.e., the eigenstates of H are close to being 
eigenstates of local operators. It is interesting to compare 
the peaks at cc > 0 which are due to local fluctuations in the 
eigenstates. Clearly, Mayer = 0 does not show any features 
because the VUMPO is diagonal in . When additional lay¬ 
ers of unitaries are taken into account, the peak structure of 
A{uj) is well approximated. The agreement in both the fre¬ 
quencies and the intensities rapidly improve with increasing 
Mayer, and the results match almost perfectly for lU = 16. 
Note that despite the extremely strong disorder, simply ap¬ 
proximating the eigenstates as product states fails to capture 
any of the interesting features. 

Comments on accuracy: We have presented some evi¬ 
dence above for the accuracy of the VUMPO obtained by 
our method. It remains to establish more precise theorems 
on what values of Mayer it would take to calculate various 
physical quantities to a specifled accuracy. In a step in that 
direction, PC have looked at the bond dimension needed to 
ensure that the smallest singular value in the Schmidt decom¬ 
position across any cut in U is less than some flxed e. This 
ensures that the discarded weight on truncating the unitary to 
bond dimension D is small. They found a slow growth of the 
^min needed to achieve a desired e with L. In the absence of 
rare resonances or Griffiths regions, T)min would presumably 
saturate at a fixed 0(1) value for a fixed error density indepen¬ 
dent of system size implying that we would need only 0(1) 
layers to represent the entire spectrum to the desired accuracy. 
As it is, with the resonances/Griffiths regions present, Omin 
is expected to grow as poly(I/) whence Mayer will grow log¬ 
arithmically. We should note however, that the PC criterion 
is not without its problems for the truncation they would em¬ 
ploy causes loss of unitarity. For example, let us return to our 
spectral function computation above but this time we first ob¬ 
tain the exact 2^x2^ dimensional unitary that diagonalizes H 
and then compress it to an MPO of a given bond dimension D. 
We do this by iteratively maximizing the “overlap” of an MPO 
with a fixed bond dimension with the “most local” diagonal¬ 
izing unitary obtained by following the PC prescription [16]. 
Because the compression scheme only minimizes the distance 
between the exact and the approximated unitary with respect 
to some operator norm, unitarity is not necessarily preserved. 
As seen in Fig. [^(labeled ED MPO), when compressing Upc 
to D = 16 (which can exactly represent our Mayer = 2 re¬ 
sults), the spectral functions A{uj) are very poorly reproduced. 
A reasonable agreement is only achieved for very large bond 
dimensions when the truncation error becomes negligible. 

Summary and discussion: We have introduced an algo¬ 
rithm to find a variational unitary MPO that approximately 
diagonalizes fully many-body localized Hamiltonians. Our 
method finds an approximation to all 2^ eigenstates of the 
Hamiltonian in a time that remarkably scales only linearly 


with system size! We have benchmarked the method by com¬ 
paring the results to exact diagonalization for small systems 
and studied the scaling of the mean variance as a function of 
system size. For a Heisenberg model in a strongly disordered 
field we find good qualitative and quantitative agreement of 
the obtained energies and spectral functions for a fixed Mayer 
and, importantly, rapid improvement with increasing Mayer- 
With this work we have provided a proof of principle that we 
can efficiently (i.e, polynomially in system size) perform a 
variational calculation that provides a complete diagonaliza¬ 
tion of fMBL systems. As the VUMPO encodes the entire 
set of eigenstates for fMBL Hamiltonians, many relevant ob¬ 
servables such as spectral functions and conductivities can be 
evaluated efficiently at zero and finite temperatures. 

A few comments are in order. First, it is intuitively clear 
that our VUMPOs should capture most of the structure of the 
eigenfunctions, or equivalently 1-bits, out to a fixed “light- 
cone” radius, set by Mayer- In terms of the dynamics this 
should allow accurate inclusion of local excitations on the 
same length scale and via the recently discussed connec¬ 
tion between the energy and size of many-body resonances 
(241 down to a related frequency scale. Indeed, this fea¬ 
ture can be effectively used to study different “slices” of 
the response function as more layers are added. For exam¬ 
ple, Figure shows that the exact solution in the case of 
W = S shows certain features at lower frequencies which 
are absent in the variational solution. Second, for a given 
VUMPO, one can construct ll25l a family of parent Hamilto¬ 
nians H = with the same eigenstates by picking 

different energy distributions for diagonal Hamiltonians in the 
“1-bit” basis, 

Going forward we can visualize many possible avenues for 
improving our method. Initially it may be possible to choose 
the same number of two-qubit gates in a different architecture 
(26l[27l to get a softer cutoff on the entanglement. More ambi¬ 
tiously we could allow for some two-qubit gates with a longer 
range and optimize over both the architecture of the unitary 
network, and the particular gates used. It is also possible to 
engineer the cost function to target a desired energy density 
via a pseudo-thermal weighting which could improve such fo¬ 
cused results for fixed resource use and also allow MBL sys¬ 
tems exhibiting mobilty edges to be treated. Of course the 
most desired improvement would be to run at Mayer > 2 
which is currently stymied by the exponential scaling of the 
cost function. As the diagrams to be contracted now start 
resembling 2D tensor-network graphs, algorithms from this 
field could presumably be used to improve the scaling of con¬ 
traction times. 

We thank Bryan Clark for useful comments on the 
manuscript. This work was supported by NSF Grant No. 
1311781 and the John Templeton Foundation (VK and SLS) 
and the Alexander von Humboldt Foundation and the German 
Science Foundation (DFG) via the Gottfried Wilhelm Leibniz 
Prize Programme at MPI-PKS (SLS). 
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Efficient evaluation of the cost functional 

In this section we discuss some details of how to efficiently 
evaluate the cost function Eq. Q using the MPO formalism. 

Due to the unitarity of U, the first term, \ 1 7 /;^), is 

simply If H is represented by a x dimensional MPO, 

the trace can by evaluated with a cost scaling as ^ Ld^x^ 
shown in Fig. (top); d is the dimension of the local Hilbert 
space on each site and is equal to 2 for the spin-1/2 operators 
considered in this work. 

The second term, {'ip^ | 1 ) ^, is somewhat more chal¬ 

lenging. We first “double” the system by taking two identical 
copies and form a tensorproduct with a state |t) (which is 
simply a product state of the “l-bits”), 

1?^') = El l^-r)|V’-r)|T). (7) 


Using the state \(p) and that we find that 


= {<P\H OF® 110). (8) 


FIG. 5. Diagrammatic representation of the tensor contractions re¬ 
quired to evaluate the terms in the cost function Eq. 0 The tensors 
represent the unitary and M the Hamiltonian. The back dots are 
delta functions (5a,6,c- 


This expectation value can again be evaluated efficiently using 
the MPO formalism as demonstrated in Fig. [^(bottom). Given 
that D > X > d, the most expensive part of the contraction 
scales as r\j LD^X^d^- 






















